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Figure 1: Alternative representations of a solid shape: (a) a surface representations, 
(b) a medial skeleton representation, (c) a curve skeleton representation. 



Abstract 

When representing a solid object there are alternatives to the use of traditional 
explicit (surface meshes) or implicit (zero crossing of implicit functions) meth- 
ods. Skeletal representations encode shape information in a mixed fashion: they 
are composed of a set of explicit primitives, yet they are able to efficiently encode 
the shape's volume as well as its topology. I will discuss, in two dimensions, how 
symmetry can be used to reduce the dimensionality of the data (from a 2D solid 
to a ID curve), and how this relates to the classical definition of skeletons by Me- 
dial Axis Transform. While the medial axis of a 2D shape is composed of a set of 
curves, in 3D it results in a set of sheets connected in a complex fashion. Because 
of this complexity, medial skeletons are difficult to use in practical applications. 
Curve skeletons address this problem by strictly requiring their geometry to be 
one dimensional, resulting in an intuitive yet powerful shape representation. In 
this report I will define both medial and curve skeletons and discuss their mutual 
relationship. I will also present several algorithms for their computation and a 
variety of scenarios where skeletons are employed, with a special focus on ge- 
ometry processing and shape analysis. 
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1 Introduction 



Three-dimensional models of solid objects are commonly used in a large variety of 
disciplines: computer graphics, medical imaging, computer aided design, visuaUza- 
tion and digital inspection, etc. The information contained in these models typically 
encodes whether a particular region of space is inside or outside the object (i.e. im- 
plicit/volumetric representations) or the interface between these two regions (i.e. ex- 
plicit/surface representations). While they are invaluable in capturing details of the 
shape at an arbitrarily fine scale, more compact and expressive representations are 
often needed. 

The advantages of compact representations are many. From a computational stand- 
point, we would like to represent the core information of our data in the most space- 
efficient way possible. This is beneficial as it allows algorithms that require exhaustive 
analysis of the input to work on much smaller datasets. The performance and quality 
of their results will depend on the conciseness and expressiveness of their description. 
On the other hand, many algorithms require some form of user intervention, like in 
editing or visualization of complex data. By providing the user with a compact and 
meaningful abstraction of the object, the complexity of this interaction can then be 
greatly reduced. 




Figure 2: Three planar solid objects Oa, Ob, Oc and their respective axes of symme- 
try. The planar symmetry relations (dotted green lines) between a few pairs of points 
(green dots) have been highlighted. 

In order to be able to generate a compact representations we can observe that most 
objects present a certain degree of redundancy in the information they convey. Indeed, 
in a geometric setting, this redundancy can be measured by the symmetry of the object. 
Intuitively, by removing the redundant information from an object that can be inferred 
through symmetry, we should be able to obtain more compact representations. 

In an object, symmetry presents itself at different levels and by exploiting it we can 
create different levels of shape abstraction. For instance, consider the symmetries of 
the 2D solid object provided in Fig. 2. For simple shapes like Oa, a single symmetry 



axis captures the global structure of the shape. However, such a global relationship is 
not always suitable. For example, by performing a simple quasi-rigid articulation of 
Oa we can obtain Ob, where a global symmetry axis is now insufficient. To address 
this problem we could use multiple, more locahzed, axes of symmetry. However, 
when the shape undergoes even more complex types of articulation, like Oc, even the 
use of a limited set of axes of symmetry becomes too restrictive. The solution to this 
problem can be found by considering symmetry at its finest level. Each pair of points 
in the shape is linked by a infinitesimal symmetry relationship and their centres of 
symmetry can be linked together to form a curvilinear axis commonly referred to as 
the skeleton of the shape. 

Historically, the extraction of skeletal representations can be dated back to Blum's 
work on the Medial Axis Transform [Blum 1967], which tracks the loci of maximally 
inscribed circles in two dimensional shapes. Although the MAT of a shape is unique, 
many alternative definitions exist, which I will discuss in Section 2.1. Each of these 
definitions highlight different properties (Section 2.2), which in turn can be exploited 
for its computation (Section 2.3). Of all these alternatives the most noteworthy in the 
context of this discussion is the one provided by [Blum 1973], which re-formulates the 
medial skeleton as the Symmetry Axis Transform of the shape: an attempt to identify 
infinitesimal (i.e. local) reflectional symmetries. 

In two dimensions, we can use the infinitesimal reflectional symmetries encoded by 
the MAT to create curvilinear skeletons. This simple network of curves provide us 
with the desired concise shape representation. Unfortunately, when we extend our 
attention to the three dimensional objects it retains its curvilinear characteristics only 
when we restrict our consideration to networks of tubular structures. For more gen- 
eral shapes, the medial skeleton, while still computable, is composed of sheets (two 
dimensional manifolds) stitched together in a complex fashion [Giblin and Kimia 
2004]. Even though medial representations of 3D objects are still useful in many 
applications (Section 2.4), there are situations in which a simpler curviUnear repre- 
sentation is beneficial. These curvilinear representations of 3D shapes are commonly 
referred to as curve skeletons and will be the subject of Section 3. 

Apart from medial skeletons, where unique and precise definitions (Section 2.1) are 
available, there are multiple ways to define curve skeletons and their computation 
(Section 3) is consequently application dependent. Typically we can think of the 
skeleton as a graph-like coding of the shape's essential structiu'e (Section 3.1). This 
encoding attempts to represent shapes in ways that agree with human intuition: by 
representing connected components and the way they connect to form a whole. Thus, 
curve skeletons can be thought of as shape descriptors: a compact and expressive 
representation suitable for solving tasks which are computationally expensive if per- 
formed on large data sets. 




Figure 3: Medial Axis Transform essentials: (a) a planar solid object (gray), its 
medial axis (white) and one of its maximal inscribed balls (red); (b) the superposition 
of a subset of medial balls and its ability to capture the original geometry; (c) the 
branch structure of the medial axis can be used to decompose the shape in meaningful 
parts. 

2 Medial skeletons 

The concept of medial skeletons originated in the context of two dimensional shape 
understanding. It was created as a way of reducing the large amount of information 
carried by a shape down to a "skeleton" of crucial information that can be more readily 
assimilated. In Figure 3 an example of medial skeleton for a planar shape is shown. 
Note how the medial skeleton is able to effectively capture the parts that compose the 
shape and the way in which they are interconnected to form the whole. 

Medial skeletons can be defined in many ways, as discussed in Section 2.1. As the 
medial skeleton of a shape is unique, these definitions are equivalent and reporting 
just one of them unequivocally defines medial skeletons. Nevertheless, these alterna- 
tive formulations are very informative: they reveal essential shape properties which I 
will further discuss in Section 2.2. These properties will offer key insights that will 
help provide the foundations for the creation of algorithms for skeleton extraction 
(Section 2.3) and their practical application (Section 2.4). 

2.1 Defining medial sl<eletons 

In this section we will formally define medial skeletons of a shape. In this context a 
shape is nothing but a solid object O C R" with boundary S — dO. Note that while I 
will define them and discuss their properties predominantly in two dimensions, these, 
unless otherwise noted, generalize to three dimensions. In we will refer to the S 
as the "contour" of the shape while in R"^ it will be simply referred to as the "surface" 
of the object. Also note that while in general the medial skeletons can be defined both 
for O and its complement space R" \ O, our focus will be on the skeletons of O: the 
internal medial skeleton. 

Maximally inscribed balls The original idea of medial skeletons was introduced 
in the two dimensional domain by [Blum 1967]. Blum proposed that the skeleton of 




Figure 4: The four possible definitions of medial axis, in lexicographical order: (a) 
only the centres of maximally inscribed balls, colored in red, are part of the medial 
axis; (b) the medial axis as the shock graph of grass-fire evolution, the location of 
the boundary at different times is illustrated by thin white lines; (c) only points on 
the medial axis have more than one corresponding point on the shape boundary; 
(d) the medial axis links boundary points with an infinitesimal reflectional symmetry 
relationship. 

an object O can be extracted by fitting maximally inscribed balls^in O and noting the 
locus of their centres (see Fig. 4-a). By augmenting the location of these centres M 
with the radius of the corresponding spheres TZ we obtain the Medial Axis Transform 
[M,TZ]^MAT{0). 



Definition 2.1. The Medial Axis Transform of O is the set of centres ^A 
and radii TZ of all the maximal inscribed balls in O. 

Note that since we only consider balls which are maximal, many inscribed balls and 
their associated loci will simply be discarded. The representation that we will obtain 
will consequently be sparse: the medial axis of O C R" has at most dimensionality 
n — 1. In addition, by superimposing medial balls of radius TZ at locations M. we can 
recreate the original shape O. This allows us to understand how the MAT is indeed a 
transform as it admits invertibility. 

While Blum's original definition provides us with much insight on its capabilities, it is 
hard to think in terms of maximally inscribed baUs. For this reason a more commonly 
known definition is the one based on the grassfire analogy. 

Grassfire analogy Imagine O c as if it were a patch of grass whose whole 
boundary S is set on fire at t = 0. The fire will propagate isotropically from the 
boundary towards the interior of the shape with uniform speed along the normals n. 
At certain locations, the fire fronts coming from different parts of the shape will meet 
and quench, thus defining a shock graph [Kimia et al. 1995]: 

Definition 2.2. The Medial Axis Transform of O with boundary S is 
given by the shock graph of the motion S = —n and the time t at which 
each shock is formed. 

This definition not only stands at the core of many medial axis computation algo- 
rithms, which we will analyze in Section 2.3, but it also allows us to intuitively under- 
stand why these skeletons are called medial. As two different sides of the boundary 
move at the same speed, their quenching location is located deep in the object, along 
its centreUne. 

Maxwell set As the grassfire propagates isotropically, quench points will always 
be equidistant from the boundary. Consequently medial points are always associated 
with at least two (Euclidean) closest points on S. This very property lies at the core 
of the Maxwell set definition of Medial Axis Transform [Mather 1983]: 

Definition 2.3. The Medial Axis Transform is the set of locations M 
intemal to the object with more than one corresponding closest boundary 
point and their distance TZ from the boundary <S. 

This definition has a key importance in the geometric setting. In Section 2.3.4 we will 
discuss how this lies at the basis of medial extraction algorithms based on the Voronoi 
diagram. 

Symmetry axis Let us momentarily restrict our attention to objects with contin- 
uous boundaries. If we take a points x S R^, then the closest point on the boundary 



s* e iS will define a ball Bx,r centred at x of radius r = \\x — s* {x) \\. This ball can 
be proved to be tangent to 5 at s* [Amenta et al. 1998b]. If we restrict our attentions 
to medial points, which have at least two closest boundary points, then the associated 
medial balls will be bi-tangent to S: 

Definition 2.4. The Medial Axis Transform of O is the set of centres 
A4 and radii 7?. of all the inscribed balls in O which are bi-tangent to its 
boundary <S. 

Interestingly, the locus of bi-tangent balls was introduced by [Giblin and Brassett 
1985] under the name of symmetry set. Consequently, the medial axis can be seen as 
a subset of the symmetry set, where only inscribed balls are to be considered. This 
simple observation underlines how skeletons can be thought of as a representation of 
symmetry axis as discussed in Section 1. Furthermore, note that while this equiv- 
alency required smooth boundaries, it extends to non-smooth ones by considering 
bi-tangency in the general sense. 

2.2 Medial axis properties 

The definitions from Section 2.1 gave us much insight on what medial skeletons, 
or better. Medial Axis Transforms are. First, we understood why it's medial, as it is 
located in the middle of the shape. Then we understood why it's an axis, by observing 
its interpretation as an infinitesimal symmetry axis and its correlation to the symmetry 
set. But most importantly we understood it is a transform as it is unique and fully 
invertible. Consequently, the Medial Axis Transform can be thought of as a dual 
representation of shape, which exhibits particular properties of the shape. In this 
section we will focus on three core properties: 1) its abiUty to compactly represent 
topology, 2) the association between medial branch and shape part and 3) the ability 
to encode local shape volume. 

In two dimensions, skeletons have received much attention as they are very effec- 
tive at capturing the topological properties of the shape. Indeed, it can be shown 
how a shape and its medial axis are homotopy equivalent [Lieutier 2003]. Although 
both volumetric and boundary representations have this ability, the skeleton is able 
to represent topology more compactly by means of a lower dimensional construct. 
Remember that for a R" solid, the medial skeletons will have at most dimensionality 
n — 1. Consequently, for two dimensional objects this will result in one-dimensional 
(curviUnear) medial skeletons. In a similar fashion, a three dimensional tubular shape 
can be represented by a medial skeleton as a collection of centrelines that split and 
join to encode the object topology. However, when we move our attention to 3D 
solids, medial skeletons wiU not be curvilinear, as they can contain up to (3 — 1)D 
construct: a collection of two-manifolds, curves and points in space (see Figure 1- 
b). For this reason the representation of topology in three dimensional data is more 
connmonly performed by curve skeletons (see Figure 1-c) discussed in Section 3. 



Figure 5: (a) Medial axis of a simple rectangular polygon, (b) Inserting a new feature 
creates localized changes in the skeletal structure, (c) The classical medial axis is 
unstable w.r.t. small scale features (i.e. noise) 



Medial branches and the "curse of instability" By observing Figure 3 we can 
see how the medial skeleton does not exclusively capture the topology of the shape. 
Indeed, there would be no topological difference in between the tuning fork and any 
other planar shape homotopic to a disk. In fact, the medial is composed of a set 
of branches which are descriptive of the shape structure - it decomposes a shape in 
parts and they way in which they are interconnected (Figure 3-c). In order to better 
understand the limitations of medial skeletons, it is important to see what causes the 
creation of skeletal branches. To analyze this behaviour we can restrict our attention 
to planar polygons, whose medial axis transform simply consists of straight lines 
and parabolic arcs. In this simplified scenario, illustrated in Figure 5-a, we have a 
straightforward relationship: a medial branch can be found for every convex vertex 
of the boundary. This is indeed a cause of concern, since as illustrated in Figure 5- 
(b,c) polygonal objects with very similar geometry can result in very different medial 
skeletons. Ideally we would desire to have a branch for any "important" geometric 
part/feature of the shape. Here, the concept of importance is to be intended in terms 
of relative geometric scale w.rt. other shape features. Note however that the classical 
medial skeleton is scale insensitive: the addition of an arbitrarily small geometric 
feature results in large, yet local, change on the medial axis. This property (or lack 
thereof) is known in the literature as the stability of the medial axis [Attali et al. 2009] 
and will be the subject of Section 2.3.5. 

Correspondence and local thickness By observing Figure 4-c we can perceive 
the existence of a very natural link between points on medial axis surfaces and points 
on the object boundary surfaces. Definition 2.3 allowed us to understand the first 
of these links: for each point on the medial axis there is one or more correspond- 
ing points on the boundary of the object. It is interesting to note that for at least 
smooth boundaries the opposite relation is quite different as for each point on 
the object boundary there is exactly one corresponding point on the medial surfaces. 
This observation is extremely important because it allows us to attach to each of the 
boundary points the size of the corresponding medial ball B. The injectivity of the 
correspondence allows us to transport the scalar function TZ defined on M onto the 
shape surface S. This scalar function defined on the surface is extremely important 



as it provides a formal way to define local shape thickness whose applications will be 
explored in Section 2.4. 

2.3 Medial axis computation 

A large body of work exists on medial axis computation. In this section I will briefly 
describe the core idea behind these solutions. In geometry processing objects are 
most commonly represented by their boundary. For this reason I wiU introduce and 
describe methods that are progressively more suitable for these representations. 

2.3.1 Exact medial computation 

Regardless of the representation, if we are given an exact description of the shape then 
it should also be possible to compute an exact medial axis representation. Typically 
these exact computations follow the intuition behind the Maxwell set definition. For 
example if the boundary of an object is represented as polygons, then we can extract 
the medial axis by computing their bisectors as shown by [Lee 1982; Held 2001] in 
two dimensions and [Culver et al. 2004] in three dimensions. Note that the extraction 
of the exact medial axis is not restricted to piecewise linear representations. Recently 
the authors of [Tzoumas 2011] demonstrated how to compute the medial axis of ra- 
tional curves in the plane, while [Musuvathy et al. 201 1] addressed a similar problem 
in the three dimensional domain. 

The extraction of exact medial representations, albeit possible, suffers from two ma- 
jor shortcomings. First of all, none of the methods above is able to deal with large 
datasets, like the ones typically found in most practical applications. This complexity 
is due to the fact that to compute the exact medial axis we have to solve systems of 
high degree algebraic equations. But even more importantly, in most scenarios geo- 
metric shapes are not known exactly. Shapes are indeed represented by different kinds 
of approximation, thus exactly extracting the medial axis is often urmecessary. For 
these reasons we will now move our attention to algorithms for approximate medial 
axis extraction and discuss the quality of this approximation. 

2.3.2 Topological thinning 

In the discussion of Definition 2. 1 I highlighted the fact that the medial axis is a 
compact structure, with a much smaller volume than the original object, and how this 
structure is necessarily contained within its boundaries. If we represent our shape as 
a set of voxels sampling its interior, then a simple idea to compute medial axis is to 
progressively reduce the volume from the outside-in until a "thin" axial representation 
of unitary thickness is obtained. This criterion lies at the base of topological thinning 
methods Uke the one in [Palagyi and Kuba 1999]. 

An atomic reduction in volume can simply be achieved by removing a voxel from 
the set. This removal has to be done carefully so not to alter the overall topology 



of the solid. Therefore, algorithms in this class typically have to specify complex 
local measures of connectedness for both interior and exterior. Also, as we generally 
desire skeletons that are medial within the shape, this removal needs to be performed 
isotropically: by simultaneously removing voxels from the whole surface. However, 
in highly discrete settings, perfect isotropy is hard to achieve. Consequently, most 
algorithms suffer of the drawback that the structure of the resulting set depends on the 
voxel processing order, often resulting in skeletons which are non smooth. In addition, 
in regions whose thickness is expressed by an even number of pixels/voxels, the set 
of centres of maximal balls is two voxels wide. In these regions further thinning is 
possible but this often results in an even rougher axis. 

For all these reasons, thinning methods are far from optimal for computing medial 
axis of detailed surface geometry as we would incur in two levels of approximation: 
in the surface to volume conversion as well as in the computation of the axis. 

2.3.3 Distance fields 

Rather than working directly on a discretized representation of the surface it is of- 
ten convenient to represent the shape by means of an Euclidean implicit function: a 
real valued function that stores the Euclidean distance from the surface to any point 
in space. As illustrated in Fig. 4-b the iso-lines of the distance field correspond to 
positions of the boundary undergoing normal motion. Exploiting this observation, 
distance field methods attempt to estimate the medial axis by analyzing the local prop- 
erties of the distance transform. Of particular interest are its ridges: the regions where 
its derivatives are multi-valued, which correspond with medial shock graphs. 

Once the distance transform is computed, the extraction voxel along its ridges can 
be performed in a number of ways. In [di Baja and Thiel 1996] the authors first 
detect critical points in the distance field and then connect them to produce medial 
skeletons. Altemative methods identify medial voxels by thresholding some scalar 
function computed from the distance field. These function generally describe the 
Ukelihood of the voxel to be located on the ridge of the distance transform, like the 
average inward flux in [Siddiqi et al. 2002] or the object angle extension of [Shah 
2005]. In many scenarios a simple thresholding is not enough as it would corrupt 
topology. For this reason, distance field methods are often coupled with the topology 
preserving thinning procedures discussed in Section 2.3.2. 

Even though implicit functions are most commonly stored on uniform discrete grids, 
they offer a much better shape approximation quality than simple voxelization, as 
typically the distance to the surface can be interpolated within each voxel. Further- 
more, while the complexity for the computation of the Euclidean distance transform 
by fast-marching [Sethian 1996] is typically 0(n log n), it is often acceptable to use 
other computationally cheaper types of distances. For example, in [di Baja and Thiel 
1996] by using city-block distance, we are able to reduce the complexity of this pre- 
processing operation down to 0{n). 



While the domain is represented more accurately, some of the limitations of thinning 
methods apply to distance transform methods as well. The medial axis is represented 
by a set of discrete voxels, and it is still possible to find regions of non-unitary thick- 
ness, making it a less than optimal choice for medial axis computation of surface 
representations. 

2.3.4 Voronoi methods 




Figure 6: (a) The Voronoi diagram of six points in the plane, (b) The Voronoi diagram 
of points sampled from the boundary of a shape. Note how Voronoi vertices and edges 
fully contained in the shape boundary form an approximation of the medial axis. 

The Maxwell set interpretation of the medial axis addressed in Definition 2.3 provides 
us with what we need to formulate the medial axis computation as a Voronoi diagram 
problem. This connection is immediately understood in two dimensions, as a Voronoi 
edge (vertex) is the locus of points that has two (three) closest neighbors. As the 
Voronoi diagram input is a point set, the boundary S of the object O typically needs 
to be sampled. Once a sampling is obtained, we can simply compute its Voronoi 
diagram V. As illustrated in Figure 6, the medial axis can then be approximated by an 
appropriate subset of V. In this section we will discuss the importance of the sampling 
and how to choose a proper subset to obtain good medial approximations. 

It was Blum in his original work [Blum 1967] that observed how, for a dense enough 
sampling of a smooth planar curve, the corresponding Voronoi diagram appeared to 
approximate A^. As illustrated in Figure 7, the denser the sampling of the curve, the 
better the medial approximation achievable. However, it was not until much later that 
[Brandt and Algazi 1992] proved that, given a dense enough uniform 5-sampling^, 
the subset of Voronoi fully contained within S correctly approximates the medial 
axis with a convergence guarantee. This result was then extended by [Attali and 
Montanvert 1997] where the requirement of an oracle knowledge of S was removed, 
and a proven topologically correct medial approximation was shown possible from 
a sampling of the surface only. Further extension was provided by [Amenta et al. 
1998a], which suggested how uniformly sampling the boundary is often unnecessary 
and an adaptive e-sampling^^can be used. 



Figure 7: The Voronoi diagram of a sampled boundary with increasing uniform sam- 
pling density. The Voronoi vertices and edges completely enclosed within the bound- 
ary approximate the medial axis. As tlie density increases, the approximation quality 
improves. Notice how a minimum sampling distance is necessary to obtain skeletons 
that are topologically equivalent to the shape. 




Figure 8: (a) A solid object in B? represented as a triangulation of its surface, (b) 
The internal subset of Voronoi vertices, (c) The internal subset of Voronoi poles, (d) 
Given a uniform surface sampling, the X-medial axis, a subset of the Voronoi facets, 
approximate the medial surfaces. 



Unfortunately, when we expand our interest to three dimensional objects, the results 
of [Attali and Montanvert 1997] or [Amenta et al. 1998a] do not directly hold. The 
problem is that, even for arbitrarily fine samplings, the Delaunay triangulation dual 
to the Voronoi diagram often presents sliver tetrahedra. These sUver tetrahedra corre- 
spond to Voronoi vertices, which neither fall in proximity of the medial axis nor are 
related to any prominent feature of the surface. To address this issue, [Amenta and 
Bern 1998] approximated the medial axis by a more restrictive subset of the Voronoi 
diagram: by only considering the Voronoi poles^ . The validity of Voronoi poles for 
medial approximation was formally verified by [Amenta et al. 2000], where the au- 
thors showed that for an e-sampled manifold, the poles approach the medial axis 
of the shape as e vanishes. The validity of Voronoi poles is intuitively illustrated in 
Figure 8 where we visually compare their approximation power to the one of internal 
Voronoi vertices. 



Figure 9: (a) A solid planar object and its medial axis, (b) Noise is added to the object 
boundary resulting in many spurious branches of the medial axis, (c) Filtering medial 
loci by object angle (X-medial axis [Chazal and Lieutier 2005b]) captures features 
across different scales but results in large topological changes, (d) Filtering medial 
loci by the distance of the corresponding surface points (^-medial axis [Amenta et al. 
2001 ]) retains topology but removes small scale features before getting rid of all the 
spurious branches, (e) The global approach of the scale axis transform is able to 
remove the noise while retaining small scale features. 

2.3.5 Stable medial axis 

When using medial representations for shape processing, we would like for similar 
shapes to have similar medial representations. However, while the Medial Axis Trans- 
form provides a perfectly invertible (i.e. dual) representation, a major shortcoming is 
that, as discussed in Section 2.2, it is not stable with regards to small geometric vari- 
ations. This issue is particularly important when we consider samples that are noisy, 
as even slight amounts of noise can produce long spurious branches in the MAT. This 
noise can be caused by a measurement process, when samples have been acquired 
by scanning the surface of real a world object, or by the fact that we are sampling a 
discretization of the shape. Here, I will review a few selected alternatives for noise 
management while referring the reader to [Attali et al. 2009] for a more extensive 
coverage of the topic. 

Typically, a simpler medial axis can be obtained by performing a smoothing opera- 
tion, however, controlling the right amount of smoothing has been shown to be a very 
challenging problem [Attali et al. 2009]. For this reason, a large majority of the lit- 
erature addresses the problem by considering the stability of a medial point either in 
terms of its local or global properties. 

Local filtering criteria One common way of creating medial skeletons is by prun- 
ing its branches according to some local criteria. A traditional choice lies in using 
some local measurement of a medial sample like the corresponding pair of points 
on the surface of Figure 4-c. For example, [Chazal and Lieutier 2004; Chazal and 



Lieutier 2005b; Chaussard et al. 2009] proposed the A-medial axis, where medial 
samples whose corresponding surface samples have a circumradius larger than A are 
discarded. For particular choices of the A parameter [Chazal and Lieutier 2005a], the 
medial axis approximation is homotopic to the shape, and its approximation quality 
is provably convergent. Unfortunately, as it is illustrated in Fig. 9-b, this criterion 
doesn't allow the capture of details across different scales, as important features of 
the shape are removed before eliminating all of the noise. An alternative choice was 
presented by [Attali and Montanvert 1996; Amenta et al. 2001; Dey and Zhao 2004; 
Sud et al. 2005], where the authors use the spoke aperture angle as a filtering cri- 
terion. However, while this approach is able to retain features at different scales, it 
might result in an axis whose topology is drastically different as we change the thresh- 
olding parameter as it is illustrated in Figure 9-c. The common problem is that since 
these solutions are highly local it allows them to be efficient and simple to analyze 
but unable to capture the multi-scale structure of the shape. 

Global filtering criterions Following a more global approach, we can compute the 
multi-scale importance of a medial branch based on volumetric approximations of 
the shape like in [Tarn and Heidrich 2003] or by analyzing its global connectivity to 
maintain topological equivalence to the shape [Sud et al. 2005]. Altematively, we 
can address the problem in a statistical sense by analyzing the resilience of a medial 
branch in a collection of shapes like in the works of [Styner et al. 2003] and [Ward 
and Hamarneh 2009]. However, because of its simphcity, the most interesting global 
solution to the problem is the Scale Axis Transform proposed in 2D by [Giesen et al. 
2009] and successfully applied to 3D data by [Miklos et al. 2010]. In the scale axis, 
the global relationship between parts is obtained by scaling the medial balls of the 
shape by a factor 5. Small medial balls, associated with small surface features and 
noisy medial branches, undergo a growth that is smaller than nearby large medial 
balls (i.e. larger nearby features). Following Definition 2.1, whenever one of these 
scaled balls is not maximally inscribed, it is removed from the set. A stable medial 
axis is then computed as the exact medial axis of the remaining set of de-scaled balls. 
As illustrated in Figure 9-e, the quaUty of the resulting medial axis is far superior to 
the one of more local methods. 

2.3.6 Approximate medial measures 

While the methods discussed up to this point provably produce convergent approx- 
imations of the medial axis, there are situations in which it is possible to compute 
coarser estimates. For example, in situations in which only a rough approximation of 
the medial surfaces is desired, the shape surface triangulation can be recycled like in 
the manifold medial approximation introduced by [Hisada et al. 2002]. Alternatively, 
in situations in which only an estimate of local thickness is required, the relatively 
expensive medial computation can be replaced by an (embarrassingly parallel) ray 
casting procedure [Shapira et al. 2008]. The distance field methods of Section 2.3.3 



can exploit coarse-to-fine spatial discretizations, where efficient yet not-as-rigorous 
algorithms for thresholding can be used [Bloomenthal and Lim 1999]. Differently 
from Section 2.3.5, the generation of stable representations can be aided by the use 
of templates offering pre-determined connectivity [Pizer et al. 2003] or by softening 
some of the strict requirements imposed by the medial axis transform [Pizer et al. 
2011]. 

2.4 Medial axis appiications 




Figure 10: (a) The approximation (left) of the local radius in [Shapira et al. 2008] 
and its application to segmentation (right), (b) Two shapes (on the main diagonal) 
are blended into each other by exploiting the volumetric structure of the medial axis 
in [Blanding et al. 2000]. (c) An object represented as a triangulated surface is 
decomposed by [Hubbard 1996] in a set of spheres to accelerate collision detection. 



The properties of the medial representation do not only lay the foundations for its 
computation but also give us the fundamental intuition to understand the many appli- 
cations it has in practical scenarios. 

The volumetric capabilities of the MAT has been exploited in a number of ways. In 
the medical domain, it has been used as a way of measuring local volumetric proper- 
ties of the shape [Naf et al. 1997], as well as for segmentation by fitting volumetric 
templates and statistical analysis [Pizer et al. 2003]. Alternatively, in 3D the local vol- 
ume can be clustered to segment the shape in patches of similar thickness as shown 
by [Shapira et al. 2008], used as a prior for hole filling [Tagliasacchi et al. 201 1] or 
to facihtate shape modeUng of organic shapes [Angelidis and Cani 2002]. The me- 
dial positioning of its curves can be exploited as an alternative to implicit functions 
to compute shape blends [Blanding et al. 2000], as well as for navigation purposes 
[Paik et al. 1998], while their interpretation as a set of inscribed spheres can be used 
to generate a simplified model for collision detection [Hubbard 1996]. In 3D, the 
boundaries of medial surfaces have been shown to be associated with important per- 
ceptual features [Hisada et al. 2002], which can be exploited to reduce the complexity 
of registration tasks [Chang et al. 2004]. As suggested by [Storti et al. 1997], the gen- 
eration of stable medial axis can be exploited in the generation of shape simplification 
schemes from a volumetric point of view like in [Tam and Heidrich 2003], while the 
continuous connection existing between medial surfaces and shape surface [Shaham 



et al. 2004] describes the volume in a way that can be exploited for the generation of 
hexahedral meshes [Tchon et al. 2003]. 

Perhaps the most important characteristic of the medial axis is that, in two dimen- 
sions, it is able to reduce the shape to a network of ID curves. For 2D shapes, this 
graph-like representation lies at the core of a large set of algorithms that address the 
important problems of shape retrieval and correspondence [Siddiqi et al. 1999]. Un- 
fortunately, these techniques do not directly apply to 3D shape because, as discussed 
in Section 2.2, their medial axis in the general setting is not composed of ID curves. 
Thus, to be able to apply these methods with minimal modifications, we need to be 
able to generate curvilinear representations for 3D shapes. These curvilinear repre- 
sentations are called curve skeletons and their computation and applications will be 
the subject of Section 3. 

We will conclude this section by describing two core applications of medial axis. 
These applications are possible thanks to the fact that the MAT offers a fully equiva- 
lent (i.e. dual) representation of the surface of an object. 

2.4.1 Surface reconstruction 



Figure 11: (a) An example of the multi-scale e-sampling; note how higher sampling 
density is required in thin regions, measured by the local feature size, (b) The me- 
dial axis is used to discard Delaunay edges connecting opposite sides of the shape's 
boundary in [Amenta and Bern 1998]. (c) An acquired point cloud (left) is recon- 
structed into a surface and its medial axis ( right) by [Amenta et al. 2001 ]. 

Given a sampling of the surface of a soUd object, the problem addressed by surface 
reconstruction is the one of obtaining a parametrized (e.g. triangulated) representation 

of the surface. 

In this context the medial axis was shown by [Amenta and Bern 1998] to be essen- 
tial in formalizing the minimal requirements of sample density. Intuitively, we would 
like a larger sampling density in regions of high curvature, so as to be able to cap- 
ture minute surface variations. At the same time, an increased number of samples is 
needed in narrow regions to properly capture the shape's topology. To address both 
these requirements, [Amenta and Bern 1998] introduced the e-sampling, where the 
necessary local sample density depends on the local feature size: the distance to the 



closest point on the medial axis. 



The contributions of medial representations do not end here. In particular, a way to 
approach the reconstruction problem is by computing the Delaunay triangulation and 
discarding a subset of its elements [Edelsbrunner and Miicke 1994]. In [Amenta et al. 
1998b], the authors extended this approach in an intuitive and powerful way: given 
an appropriate e-sample, they show that, by discarding elements whose circumsphere 
contains a portion of the medial axis, results in a topologically correct reconstruction 
of O. This method and its essential relationship with the medial axis created the basis 
for much work in the field of combinatorial reconstruction [Cazals and Giesen 2006]. 

2.4.2 Shape manipulation 




Figure 12: (a) An example of a model produced by the medial axis based framework 
of [Angelidis and Cani 2002]. (b) A graph approximating the medial axis is used 
as an embedding domain for a template skeleton, allowing to re-use animation rigs 
and produce new poses [Baran and Popovic 2007]. (c) An example of a deformation 
achieved by the volume preserving framework of [Yoshizawa et al. 2007]. 

Developing geometry-manipulation tools for end users is a core problem in geometry 
processing. The medial axis offers an alternative to the more commonly used surface 
based methods by exploiting the volumetric information of the shape captured in the 
axis. 

In [Storti et al. 1997] the authors propose to use the medial domain as an alternative 
to surfaces for shape parameterization, thus creating the foundations for the widely 
used m-reps of [Pizer et al. 2000] and the shape modeling frameworks of [Angelidis 
and Cani 2002]. Intuitively, by sparsely sampling the medial axis, the shape can be 
manipulated by modifying the position and radius of medial loci and then reconstruct- 
ing/deforming the object surface. 

Following a similar paradigm, [Bloomenthal and Lim 1999] investigated a way of 
connecting the use of geometric (medial) skeletons to traditional kinematic skeletons 
for shape animation. Rather than directly expressing surface vertex positions as a 
function of the kinematic skeleton, the author proposes to first animate the medial 
skeleton and then reconstruct the surface from its medial representation. Further con- 
nections between these two structures were also investigated by [Baran and Popovic 



2007], where instead of attempting to model surface animation by means of medial 
axis, the authors consider the medial domain an appropriate space for fitting a tem- 
plate of a kinematic skeleton. Traditionally, the transformation of animated surface 
vertices is obtained by a linear combination of transformations defined on kinematic 
skeletons. While these weights are typically set by artists by a time consuming trial 
and error process, [Bloomenthal 2002] uses the medial axis and the radius function 
defined thereon to create a convolution scheme that is able to completely automati- 
cally create blending weights, resulting in natural-looking surface animations. 

Rather than directly modifying the medial axis, or use a kinematic skeleton as a con- 
trol structure, in [Yoshizawa et al. 2003; Yoshizawa et al. 2007], the authors investi- 
gated the use of freeform deformation of the medial siu'faces for shape deformation. 
Their solution presents two strong advantages: first of all, thanks to the volumetric na- 
ture of medial representation, more natural large scale deformations can be achieved 
as local thickness is preserved; furthermore, self-intersections, which often arise when 
large deformations take place, can efficiently be corrected by exploiting the volumet- 
ric nature of medial representations. 



3 Curve Skeletons 




Figure 13: (a) A three dimensional object; the colour coding expresses the corre- 
spondence between portions of the surface and vertices on the curve skeleton, (b) An 
example of a curve skeletal representation extracted by [Au et al. 2008]. 

Curve skeletons were introduced as a way of providing a curvilinear (a one dimen- 
sional construct) representation of three dimensional shapes. As I have thoroughly 
discussed in previous sections, this necessity came from the fact that classical medial 
axis generally produces two-manifold elements when appUed to three dimensional 
shapes. It is important to note that while medial skeletons are precisely and uniquely 
defined, the same doesn't apply to curve skeletons. There is no universally accepted 
rigorous definition of what a curve skeleton is. A notable exception is the medial 
geodesic skeleton recently defined by [Dey and Sun 2006]. However, since it is com- 
puted as a subset of the medial axis, it also inherits all of its stability issues and 
has consequently received little attention in appUcative scenarios. Without restrict- 
ing ourselves to this formal definition, we can compute curve skeletons in multiple 



ways; these different methods will be the subject of Section 3.2. The lack of a com- 
mon formal background makes the task of comparing these different algorithms quite 
difficult. An attempt to compare these techniques was performed by [Cornea et al. 
2007], where the authors defined a set of desired curve skeleton properties and then 
discussed several algorithms accordingly. In Section 3.1 we will briefly discuss these 
properties, which will find a direct use in creation of curve skeleton extraction (Sec- 
tion 3.2) and their applications (Section 3.3). Note that, similarly to Section 2, 1 will 
focus my attention to methods which are suitable for surface-based object representa- 
tion of 3D geometry. 

3.1 Properties 

The only strictly necessary property that a curve skeleton should have is that it needs 
to be thin - a one-dimensional entity embedded in the three dimensional domain. 
As in medial skeletons, it is also often required that this network of curves is homo- 
topic to the shape - that they retain the ability to accurately encode the topological 
structure of the object they represent. It is also important that curve-skeletons offer 
part-awareness: its branches should provide a one-to-one association with "impor- 
tant" geometric features on the shape. As the graph connecting branches establishes 
a relationship network, by computing skeletons at different resolutions we obtain a 
way to represent shape in a hierarchical fashion. Additionally, we would like their 
computation to be robust, in the sense that the presence of surface noise should not 
result in the generation of spurious branches. The smoothness of curve skeletons is 
also important, not only for its visual appeal, but also because it improves compact- 
ness, as smooth structures are generally representable with a smaller set of geometric 
elements. 

Furthermore, as curve skeletons are often employed as representations of solid ob- 
jects, it is often desirable for its curves to be contained within the object. In some 
scenarios, an ever more restrictive requirement on its positioning is necessary: the 
curves are required to be centred within the object, in a similar way in which MAT 
skeletons are medial within the shape. The positioning of the medial axis is also 
critical in situations that require reliability: where we ask that every point on the sur- 
face is visible from at least one curve skeleton location. It is important to note that, 
from a geometrical point of view, curve skeletons do not truly provide an alternative 
shape representation, as we cannot perfectly reconstruct the original shape given a 
curve skeleton. Curve skeletons are instead to be interpreted as shape descriptors: 
a compact way of capturing essential geometrical and topological information. In 
order to be able to transfer information from the surface to the curve skeleton it is 
necessary that proper correspondences exist between points on the surface and points 
on the skeleton. For example, these are essential for the generation of an approxi- 
mated reconstruction of the object. This can be achieved by superimposing a set of 
simple primitives (e.g. spheres, ellipsoids or by sweeping planar curves) along the 
curve skeleton. The parameters of these primitives can be inferred by analyzing the 



corresponding surface geometry. 



It is important to note that not all properties are essential in every application scenario. 
Furthermore, strictly requiring a skeleton to offer one of these properties might make 
the satisfaction of others unfeasible. As an example, it is hard to obtain centred (in a 
medial sense) yet robust skeletons, because the medial axis is sensitive to small sur- 
face perturbations. Similarly, we could achieve better reconstruction by introducing 
extra branches in the representation, albeit at the cost of losing part-awareness. 

3.2 Curve skeleton computation or "skeletonization" 

In the literature, curve skeletons are computed or "extracted" following few altema- 
tive philosophies. Some algorithms, discussed in Section 3.2.1, focus on reducing the 
medial skeleton surfaces down to one dimensional structures. In Section 3.2.2 I will 
discuss those that attempt to modify the notion of distance in order to create curves, as 
opposed to surfaces, which lie medially within the shape. In Section 3.2.3 1 will move 
my focus to solutions based on topological analysis, where the core aim is to gen- 
erate homotopy preserving skeletal representations. Other very successful methods, 
described in Section 3.2.4, focus on the fact that we are indeed looking for a curve. 
I will conclude with Section 3.2.5, where I will describe very-application oriented 
techniques, that construct curve skeletons by grouping surface elements according to 
particular properties, and then define the skeleton as the graph describing how these 
components are glued together to form the whole. 

3.2.1 Skeletonization by medial axis processing 




Figure 14: An illustration of the fundamentals of the described in [Dey and Sun 
2006]. (a) Three points ( red) on the medial axis of a cuboid; their corresponding sur- 
face points (blue), and the shortest paths on the surface between these points (black), 
(b) The medial geodesic function on the medial axis of shape, (c) High divergence 
points of the medial geodesic function define the curve skeleton. Note how only y and 
z are on the skeleton as they are associated with more than one shortest path, (d) The 
curve skeleton of a noisy hand model is also noisy as it is defined as a subset of the 
medial axis. 

The medial axis of a solid object conveys much of the essential structure of the shape. 
However, the medial skeleton is geometrically complex as it is composed of a set of 
surfaces. A simple approach to compute a curviUnear representation of the shape is to 



prune the medial axis surfaces down to a linear structure. Many heuristics can be used 
to prune the medial sheets to a line-form; for example, the authors of [Lee et al. 1994] 
use a thinning method similar to the one described in Section 2.3.2. A very interesting 
solution along these lines is the one proposed by [Dey and Sun 2006], where the 
authors proposed to extract the skeleton as the shock graph of the Medial Geodesic 
Function, a scalar function defined on the medial axis surfaces. Given a point m E Ai 
and its two corresponding surface points si , S2, this function is defined as MGF(m) = 
ds{si,S2), where ds{si,S2) is the geodesic distance on S between the two points. 
The theoretical argument on the existence of this shock graph provided by [Dey and 
Sun 2006] make medial geodesic skeletons the only mathematically precise version of 
curve skeletons. Nevertheless, curve skeletons extracted in this fashion suffer severe 
limitations. Its computational complexity is large, as it necessitates geodesic distances 
between every pair of surface vertices. Furthermore, the sensitivity of medial axis to 
surface noise, discussed in Section 2.2, result in the extraction of noisy skeletons Uke 
the illustrated in Figure 14-d. 

3.2.2 Skeletonization by generalized field analysis 

There exist alternatives to processes that erode the medial axis in order to generate 
curve skeletons that are well-centred within the shape. These are based on more com- 
plex (in particular less localized) ways of measuring the distance from the boundary 
[Hassouna and Farag 2008]. The value of the field at a position in space is computed 
as the average of potential fields of many boundary samples, and the skeleton extrac- 
tion is performed by tracing curves seeded at critical points along directions of high 
divergence. Because of the inherent averaging, these skeletons are not only medial 
but often very smooth. A major problem with these methods is that they are not par- 
ticularly suitable for surface based geometry: the field needs to be represented on a 
discretized volume and, most importantly, they tend to suffer from numerical instabil- 
ities. For these reasons they will not be further discussed and I simply refer the reader 
to [Cornea et al. 2007; Hassouna and Farag 2008] for more details. 

3.2.3 Skeletonization by topological analysis 

There are situations in which the core information about an object can be captured by 

the topological structure of its skeleton (i.e. loops and joints) rather than by its geo- 
metric attributes (i.e. curve locations and local radius). In Morse theory [Shinagawa 
et al. 1991], a real function / defined on the surface <S of a shape is able to capture 
its essential topological information: the critical points of / encode the topological 
changes in the shape. More intuitively, the topological information can be extracted 
by looking at the iso-contours of / and their distribution on S [Biasotti et al. 2000b; 
NataU et al. 2011] as illustrated in Figure 15. 

In certain situations, like Figure 15-b the simple y coordinate of boundary points 
is sufficient. However, there are many ways of creating / and this choice is made 



Figure 15: (a) A geodesic distance field sourced at the marked point is created onto 

the surface and its iso-contours are shown. The barycenter of surface portions be- 
tween two successive iso-contours is used to build a curve skeleton like embedding 
of the topological graph [Lazarus and Verroust 1997]. (b) The y coordinate is used 
as the field for topological analysis. By recursively refining its iso-contours, a multi- 
resolution topological graph is created [Hilaga et al. 2001]. (c) The iso-contours of 
the average geodesic function from [Hilaga et al. 2001 ] naturally capture the extrem- 
ities of the articulated shape of a frog. 

according to the type of information we want to extract from the shape. While for a 
broad comparative analysis we refer the reader to [Biasotti et al. 2003], for skeleton 
extraction, a simple solution was introduced by [Lazarus and Verroust 1999]. Their 
idea is illustrated in Fig. 15-a, where / is computed as the geodesic distance (or one 
of its approximations) from a selected source point p. Even though the seed point can 
be selected by heuristics, it is important to note that a single seed generally results in 
the selection of a preferred slicing direction. An example of the consequences of this 
choice is illustrated in Fig. 15-a, where we can easily discern a loss of structural detail 
in the skeleton. To solve these issues, [Mortara et al. 2003] proposed the use of a set of 
seed points on protrusions, which can be identified by searching for locations of high 
large gaussian curvature, while other methods like [Hilaga et al. 2001] compute the 
integral of the geodesic distance to all other surface points (See Figure 15-b,c).The 
limitation of all these methods is the necessity to pick an heuristic, as, according 
to the choice made, unforeseen topological loops can appear in locations where the 
wavefronts collide. 

It is important to note that this class of algorithms typically offers limited geometric 
interpretation. The spatial location of a skeletal vertex (i.e. its embedding) is gen- 
erally computed a posteriori - as the centroid of its corresponding surface points. 
Consequently, these methods generally offer lower quality geometric skeletons, and 
are more focused toward providing a correct topological representation. It is very 
important to note that in appUcations where topology is the core shape feature, very 
simple (and efficient to compute) choices for / can be made, like the y coordinate of 
Figure 15-b, and these are ensured to produce homotopy preserving curve- skeletons. 



3.2.4 Skeletonization by contraction 



The fundamental requirement of curve skeletons to be a line-like representation can 
be used for their computation. A line representation can indeed be thought of as a 
zero-volume solid. Thus, in order to be able to extract a skeleton, we can let the shape 
evolve in such a way that, at each iteration, its total volume is strictly decreased. In 
order to avoid loss of detail, this volume reduction, commonly referred to as con- 
traction, needs to be performed in a controlled fashion. The method used to control 
contraction is the fundamental ingredient of the algorithms described in this section. 

Figure 16: An illustration of the volumetric contraction proposed by [Wang and Lee 
2008]. (a) The surface is first converted into a volumetric representation, (b) The 
edges of the volume are isotropically contracted, but by using the surface attraction 
constraints (blue) we can avoid the shape from collapsing to a point, (c) By using 
these attractors the collapsed volume retains the shape features. 

Volumetric contraction In [Wang and Lee 2008], the authors compute the curve 
skeleton of a volumetric model by applying a contraction procedure. If we think of 
the shape's volume as a set of voxel centres connected by edges, it is intuitive to 
think that we could reduce the shape's volume by uniformly scaling down the edge's 
lengths. Clearly this isotropic scaling would quickly and meaninglessly shrink the 
shape down to a zero- volume point as illustrated in Figure 16-b. In order to prevent 
this from happening, the authors solve a constrained problem where the locations of 
voxels originally positioned nearby the surface are required not to move too far from 
their initial location. A fundamental problem in applying this technique to surface 
models is that it requires a volumetric contraction. Nevertheless the contribution of 
[Wang and Lee 2008] is important as it built the groundwork for the very successful 
surface based contraction method. 

Surface contraction In [Au et al. 2008] the authors address the problem of vol- 
umetric contraction directly on a meshed surface representation. In the proposed 
solution they iteratively reduce the volume of the shape (i.e. contraction) by solving a 
constrained mean curvature flow (MCF). Motion by mean curvature is known to result 
in a volume loss at each iterations [Ruuth and Wetton 2003]. While performing MCF 
on the surface doesn't cause a contraction as isotropic as [Wang and Lee 2008], sur- 
face attraction constrains need to be inserted in the system to make sure that branches 
do not undergo excessive shrinkage, thus, resulting in the loss of important surface 




Figure 17: An illustration of the surface contraction proposed by [Au et al. 2008]. 
(a,b,c,d) From left to right the shape is contracted by applying curvature flow until the 
volume vanishes. To prevent over-contraction surface attraction constraints similar 
to [Wang and Lee 2008] are applied, (e) The contracted surface retains the original 
manifold topology and is collapsed to linear graph [Li et al. 2001 ]. 

features. The foundations for the success of this technique is that both the volume 
contraction motion and the feature preserving constraints can efficiently be encoded 
by a system of linear equations for which efficient solvers are available. Furthermore, 
as far as a proper Laplacian operator is used, the contraction procedure is not only in- 
sensitive to badly discretized surfaces (e.g. low quahty triangulation), but since MCF 
is essentially a smoothing operation, it is also extremely robust to noise. 

Topological contraction While the methods above tackle the problem of contrac- 
tion in a geometrical sense, it is interesting to know that it's possible to reduce the 
overall shape volume by only using topological operations. In [Li et al. 2001], the 
authors propose to use iterative edge collapse to reduce a meshed surface to a set 
of simple segments, which have no incident faces (i.e. wingless). As the resulting 
set is a simple network of edges, this structure is ensured to be one dimensional. 
Furthermore, by choosing to collapse the shortest edge first (resulting in an overall 
0{n log n) complexity), the authors experimentally demonstrate how this solution re- 
sults in appro ximatively-centred skeletons. While no strong theoretical foundations 
exist for topological contraction, this algorithm lies at the core of many skeletoniza- 
tion algorithms, as it is often employed as a post-processing step: to reduce a quasi- 
thin structure to a set of truly one dimensional primitives [Wang and Lee 2008; Au 
et al. 2008]. It is important to note that in situations in which the quasi-thin struc- 
ture is already centred within the shape, the topological contraction does not alter the 
embedding, thus, preserving centeredness. 

3.2.5 Skeletonization by property grouping 

If we are given the curve skeleton of a shape, we are provided with a way of decom- 
posing it in fundamental parts and their interconnections. Note that this relation can 
be considered in the other direction as well. If our apphcation domain allows us to 
precisely identify what constitutes a part, then we can construct curve skeletons as 
the network of curves describing their spatial configuration. There are many ways to 
define appropriate grouping criteria for skeleton generation and, in this section, I will 



Figure 18: (a) A convex decomposition allows to recycle PCA axes as skeletal 
branches [Lien et at. 2006]; the finer the decomposition, the more accurate the skele- 
tonization, (b) Local rotational symmetry of a group of samples can be measured by 
observing the arrangement of the surface normals on the gaussian image [Tagliasac- 
chi et al. 2011]. (c) The articulation of a shape is essentially described by a set of 
locally rigid transformations [Anguelov et al. 2004]; these transformations can be 
identified and exploited to create a kinematic skeleton. 

describe a few important examples. 

Rotational symmetry In Section 1 1 discussed how one way of interpreting a skele- 
tal axis is by considering it as an axis of symmetry. For two dimensional shapes, it was 
refiectional symmetry that allowed us to define medial skeletons - by grouping pairs 
of locally symmetric points [Zhu 1999]. In three dimensions a different symmetry re- 
lationship can be exploited: rotational symmetry. Nevertheless, like in Section 1, the 
extraction of global rotational symmetry [Sun and Sherrah 1997] would only allow us 
to define a single global axis. Consequently, the generation of curvilinear representa- 
tions can be achieved by considering symmetry at a local scale [Li et al. 2001]. The 
identification of local rotational symmetry (see Figure 18-b) lies at the core of a class 
of algorithms that extract skeletons by exploiting this geometric observation [Chuang 
et al. 2004; Tagliasacchi et al. 2011]. 

Segmentation Aside from symmetry, there are many ways to partition an object 
into geometrically meaningful components [Shamir 2008]. Segmentation schemes 
which are meaningful in skeletonization usually perform the shape decomposition by 
following some convexity criterion - by cutting the surface along concave regions. 
Decomposing a shape into convex components creates a segmentation that is mean- 
ingful for skeletonization, as the principal axis of a convex component (in a PCA 
sense) is a linear element that is centred and reliable [Lien et al. 2006] (also see Fig- 
ure 18-a). Another way of interrelating segmentation to skeletonization is found by 
considering the animation of articulated objects. For example, the concavity formed 
by a bent elbow gives us strong cues as to where the articulation takes place. The 
identification of these concavities is used to decompose the shape into components 
and allows the generation of hierarchical animation skeletons [Katz and Tal 2003]. 



Articulation analysis The skeleton computation algorithms considered up to now 
mostly addressed static scenarios, where the shape doesn't change in time. However 

dynamic data can also be exploited in the skeleton extraction process. In articulated 
shape animation, this dynamic data consists of several frames of an object in differ- 
ent poses. These frames can be analyzed to discover the underlying structure. The 
data contains large connected areas undergoing a similar quasi-rigid transformation. 
Following this intuition, [Anguelov et al. 2004; Schaefer and Yuksel 2007; De Aguiar 
et al. 2008] construct curve skeletons by first clustering parts of the shape undergo- 
ing similar transformations, associating a bone to each of them and connecting them 
together to create an animation rig resembling the ones used in animation (see Sec- 
tion 3.3.1). 

3.3 Applications 

Because of their simple yet informative structure, curve skeletons have found uses 
in many areas including scientific visualization, computer animation, shape analysis, 
and computer assisted medical diagnosis, etc. The coverage provided by this section 
is intentionally not exhaustive and it will be strongly biased toward applications in 
computer graphics and geometry processing. 

3.3.1 Computer animation 




Figure 19: (a) An example of a kinematic skeleton (also called animation rig) from 
MAYA is used to articulate the shape from its rest-pose to a desired one. (b) Each bone 
(i.e. edge) of the curve skeleton is associated with skinning weights that describe 
how its movement affects the surface vertices [Jacobson et al. 2011]. (c) Given a 
set of example poses (top), it's possible to automatically estimate skinning weights 
and skeletons [Schaefer and Yuksel 2007]; (bottom) this information is sufficient to 
extrapolate new poses. 

In computer animation, posing digital characters is an essential task and there are mul- 
tiple ways to achieve it. For example, we can directly manipulate the surface [Lipman 
et al. 2005] or use control cages [Lipman et al. 2008]. Arguably the most commonly 
used method is Skeletal Subspace Deformation [Lewis et al. 2000], where curve skele- 
tons mimic the way real-world bones create deformations in our own bodies. Not only 
is skeletal deformation foimd in most 3D modeling packages but its efficient compu- 



tation is so crucial that it is commonly hardware-accelerated on GPUs. In SSD, the 
animator specifies a kinematic skeleton for a rest-pose of the shape - Uke the T-pose 
of Figure 19-(a,left). To obtain a desired pose, a set of transformations (typically 
rigid) are applied in hierarchical fashion to the rest-pose skeleton. Finally, the posed 
shape is obtained by transferring the bone transformations onto surface vertices as 
shown in Figure 19-(a,right); this is achieved by performing a weighted combination 
of transformations. These weights are typically called skinning weights: a scalar field 
that encodes how much a surface region should be influenced by a particular bone 
(See Figure 19-b). 

As the generation and skinning of kinematic skeletons is a time consuming task, a 
number of algorithms has addressed the problem of automating the process. For 
example, in [Wade and Parent 2000; DeUas et al. 2007] the authors considered the 
problem of fitting templates of skeletons in digital human models. In [Baran and 
Popovic 2007], the authors generalized the approach, not only considering fitting of 
more complex skeletons (i.e. not human), but also in the more general sense of trans- 
ferring animations between different shapes. If a sequence of pre-animated meshes is 
given, a skeleton can be extracted with the techniques discussed in Section 3.2.5 and 
their skinning weights can be computed in a completely automated fashion [Schaefer 
and Yuksel 2007] as illustrated in Figure 19-c. This possibility of converting anima- 
tions to the SSD framework has been demonstrated to be very effective for animation 
compression [De Aguiar et al. 2008]. In addition, note that the generation of kine- 
matic skeletons and weights is not restricted to surface models, but elegantly extends 
to volumetric data [Gagvani and Silver 2001; Theobalt et al. 2004]. 

The applications of curve skeletons in animation are not limited to generating shape 
poses and parameters. For example, in [Lazarus and Verroust 1997], the authors 
employ skeletons as a parameterization domain suitable for shape blending of sim- 
ple shapes. Furthermore, the performance of collision detection algorithms based on 
bounding volume hierarchies can be enhanced by decomposing shapes using skele- 
tons [Lietal. 2001]. 

3.3.2 Geometry processing 

Curve skeletons have broad appUcations in the modeling and analysis of geometry. 
Their branches are part-aware and associated with important shape features, allowing 
the development of efficient segmentation algorithms for simple human-like objects 
[Xiao et al. 2003; Werghi et al. 2006], as well as for shapes with more complex struc- 
ture [Reniers and Telea 2007; Tiemy et al. 2007; Au et al. 2008]. Their ability to 
compactly represent topology has been exploited to automatically correct topological 
problems in isosurfaces of implicit functions [Wood et al. 2004] and low quality sur- 
face models [Zhou et al. 2007]. Its approximate reconstruction properties have been 
exploited for the modeling of organic shapes [Angelidis et al. 2002; Ji et al. 2010] as 
weU as for the correction of geometry in acquired data [TagUasacchi et al. 2009; Cao 



Figure 20: (a) The surface to skeleton correspondences allow to define a radius 
fiinction (top) onto the surface which can be used to segment the shape (bottom) [Au 
et al. 2008]. (b) The compact topological representation of the curve skeleton is used 
to drive a topology fixing process [Zhou et al. 2007]. (c) Multi-dimensional scaling 
can be efficiently performed on the skeleton, resulting in stretched out limbs which 
facilitate the task of left/right disambiguation commonly found in articulated shape 
correspondence [Au et al. 2010]. 

et al. 2010; Zheng et al. 2010]. Curve skeletons are also very useful to address the 
correspondence problem particularly for datasets consisting of shapes characterized 
by many protrusions. Their compactness allows to greatly reduce the input size of 
what is essentially a combinatorial problem. 

Shape correspondence Shape correspondence addresses the problem of identi- 
fying meaningful relationships between elements of two shapes. A large variety of 
methods solve the problem by directly analyzing the shape geometry [Tangelder and 
Veltkamp 2008; van Kaick et al. 2010]. On the other hand, curve skeletons compactly 
represent a shape as the graph of its components, their relationship and, in some sce- 
narios, a rough description of their geometry [Biasotti et al. 2000a]. Consequently, an 
alternative approach to establish a correspondence between shapes can be obtained by 
comparing their graphs [Hilaga et al. 2001]. This is advantageous because the size of 
graph is usually much smaller than the size of the surface (i.e. number of triangles) or 
volumetric representations (i.e. number of voxels), allowing more computationally- 
expensive algorithms to run in reasonable time. In addition, when shapes exhibit 
large variation in poses and surface details, a graph representation better captures the 
shape's overall structure than local geometric measures (e.g. curvature). 

A possible way to compare two graphs is to identify the maximal common subgraph, 
but this problem is NP-hard. Alternative methods employ graph edit distances - by 
measuring the cost of transforming one graph into the other by simple editing oper- 
ations. Although this is known to be an NP-hard problem [Zhang et al. 1996], there 
exist heuristic solutions (e.g. the graph often needs to be converted into an acyclic 
one) based on graph editing operations both in 2D (Sebastian et al. 2004] and 3D 
[Biasotti et al. 2006]. A different way to approach the problem is by noting that the 
spectral decomposition of the graph's adjacency matrix contains much information 



about the topology of a graph [Siddiqi et al. 1999; Sundar et al. 2003]. As the eigen- 
decomposition typically captures global topological information about the graph, al- 
gorithms in this class typically proceed in a coarse to fine manner. For this reason 
the skeletal graph needs to be not only cycle free but also capable of expressing the 
natural hierarchical structure of the shape. Other more direct solutions compute cor- 
respondences by directly comparing sets of skeletal nodes in either a hierarchical way 
by means of a set of heuristics [Hilaga et al. 2001], or by a carefully designed voting 
scheme [Au et al. 2010]. 

3.3.3 Medical shape analysis 




Figure 21: Examples of the use of curve skeletons in medical applications, (a) The 

skeletal centreline of the colon is extracted to allow virtual colonoscopy [Wan et al. 
2001] or its unwinding [Silver and Gagvani 2002]. (b) The skeleton of the vessel 
structure is used to perform planar reformation to visualize its untangled structure 
[Kanitsar et al. 2003]. (c) The centreline automatically generates flight paths for 
virtual bronchoscopy [Perchet et al. 2004]. 

In the medical domain, curve skeletons are well suited to describe anatomical struc- 
tures that have characteristic tubular shapes like vessels, nerves and elongated mus- 
cles [Nystrom et al. 2001; Fridman et al. 2004]. The compactness of the represen- 
tation allows to efficiently perform computational tasks like registration of partially 
overlapped vessel images [Aylward et al. 2003] and the flattening of their three di- 
mensional structure to a plane [Kanitsar et al. 2003]. Furthermore, the encoded local 
volumetric information can help in the detection of abnormalities in vascular struc- 
tures, Uke stenosis [Sorantin et al. 2002] and aneurisms [Straka et al. 2004]. Skeletons 
satisfying centeredness and reliability are also particularly important in medical ap- 
plications. Indeed, in virtual navigation [Wan et al. 2001], it is essential that the 
physician is able to examine the full internal surface of an organ [He et al. 2001]. 
In this context, appropriately extracted curve skeletons can be used to generate cam- 
era fly-through paths for the inspection of the intestine [Hong et al. 1997; Wan et al. 
2002], lungs [Perchet et al. 2004] and blood vessels [Bartz et al. 1999]. In addition, 
when the structure is so tangled that even a virtual navigation would be uninforma- 
tive, the skeleton can be used to straighten the organ to a linear structure which can 
be more easily analyzed [Silver and Gagvani 2002]. 



4 Conclusions 



In this report, I presented a detailed overview of skeletal representations as an alter- 
native (or addition) to boundary or volumetric representations for shape processing. 

First, I introduced the idea of medial skeletons as a way of capturing redundancy 
and obtaining local axial representations. I discussed different possible interpreta- 
tions, several classes of algorithms for its computation, its limitations in terms of 
stability, along with several ways to tackle this problem. For surface-based repre- 
sentations, I justified how Voronoi-based medial axis extraction methods offer clear 
advantages with regards to their competitors: not only do they offer theoretical con- 
vergence guarantees under appropriate sampling, but efficient and robust libraries for 
three dimensional voronoi computation are freely available as well. Although, only 
for two dimensional shapes, medial skeletons result in curvilinear structures, 1 illus- 
trated how the medial surfaces associated with three dimensional solids contain useful 
information exploitable in many appUcations. 

The necessity for a simpler curvilinear structure for three dimensional shapes stirred 
the discussion toward curve skeletons. While there is no unique definition for curve 
skeletons, 1 concisely described the properties that these curvilinear representations 
are commonly required to have. 1 then presented several classes of algorithms for 
skeletonization of solid shapes, focusing on those algorithms that operate on bound- 
ary representations. I highlighted how contraction methods currently represent the 
state of the art for skeletonization, because they can efficiently extract skeletons from 
low quality and possibly noisy surface representations without needing any domain- 
specific information. 1 concluded my discussion by describing a number of appli- 
cations of curve skeletons and highlighting their abiUties in describing shape from 
various applicative fields. 
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Appendix 



Maximal inscribed ball Let O be a closed connected set in R". A closed ball 
B C R" is called maximal inscribed ball in O if B C O and there does not exist 
another ball B' ^ B such that B c B' c S. 

Voronoi diagram Given a finite set of points S in R'' , for each point p in R*^ there 
is at least one point in S closest to p. a point p may be equally close to two or more 
points in S. For each point in S its Voronoi cell is defined as the subset of -R*^ of points 
closest to it than to any other point in S. The union of Voronoi cells of all points in S 
is called the Voronoi Diagram of the set S. For n = \S\ and /c = 3 the complexity is 
O(n^), but when points are well distributed on a smooth surface the complexity has 
been shown to reduce to 0{n log n) [Attali et al. 2003]. 

Voronoi pole Given a finite set of points S in R'^ and its Voronoi diagram, each 
sample point p is associated with a convex Voronoi polyhedron. The vertices on the 
polyhedron on the two sides of the boundary which are further from p are the Voronoi 
poles of p. Note that there are at most 2n poles for a sample of size n. 

Local feature size The local feature size LFS(a:) is the minimum Euclidean dis- 
tance from any point x to any point of the medial axis. 

Delta sampling (i.e. uniform sampling) A sample P of <S is a (5-sample if no point 
p e «S is further than 5 from P. 

Epsilon sampling (i.e. adaptive sampling) A sample P of <S is an e-sample if no 
point p G <S is farther than r ■ LFS{p) from M. 



